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Abstract 

Point defect migration is considered as a mechanism for aging in ferroelectrics. Numerical results 
are given for the coupled problems of point defect migration and electrostatic energy relaxation 
in a 2D domain configuration. The peak values of the clamping pressure at domain walls are in 
the range of 10^ Pa, which corresponds to macroscopically observed coercive stresses in perovskite 
ferroelectrics. The effect is compared to mechanisms involving orientational reordering of defect 
dipoles in the bulk of domains. Domain clamping is significantly stronger in the drift mechanism 
than in the orientational picture for the same material parameters. 
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I. INTRODUCTION 



Ferroelectric materials underlie restrictions in technological applications because of sev- 
eral degradation phenomena. One of these phenomena is aging which is defined as the 

Badual change of material properties with time t under static external boundary conditions 
|8|. Experimentally, the dielectric constant decreases and hysteresis loops 
alter their shape and amplitude in remanent polarization and coercive field. The dielectric 
constant decreases as the logarithm of time in an intermediate time range saturating for 
long times . The change of material properties is caused by a decreasing domain wall mo- 
bility stabilizing in an aged domain structure |lo|. In order to describe_the experimentally 
observed shifted or deformed ferroelectric hysteresis loops after aging 



the internal field 

E* has been defined as a quantity describing the strength of domain stabilization j^. Sev- 

le stabilization process, space charge 



eral mechanisms have been considered to partake in t 

formation 3, 4], domain splitting [lo|, ionic drift 11, 12, 13 1, or the reorientation of defect 
dipoles iJ] . Except for domain splitting, all mechanisms directly involve some reordering of 
point defects. Within a microstructure, three relevant locations can be identified for charge 
carrier rearrangement: the domain bulk, grain boundaries, or domain walls. For the bulk 
effect, defect dipoles reorient with respect to the direction of the spontaneous polarization 
under an electrical field or strain. For the grain boundary effect, charged point defects move 
under electrical fields originating from polarization discontinuities at the grain boundaries or 
the outer perimeter of the sample in order to compensate the fields. The same process can 
occur at charged domain walls and then becomes a domain wall effect Is], [l^. Local space 



charge is another electric driving 



17 



181 ] as well as perovskites 
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brce for ionic currents observed in LiNbOa type crystals 



20 



2l|. 



Elastic fields can provide a second driving force for defects inside domains but will not 
be treated here. For not too rigid non-charged domain walls, the localization of free charge 



carriers at a domain wa^ 



planar walls 22, 



23|, 



1 is a further possible effect entailing the wrinkling of the initially 



25(1 . Even though this is a physically interesting mechanism, it will 



not be taken into consideration here. 

Oxygen vacancies are a well known and frequent defect in the perovskite structure. They 
have been considered to play an important role in aging of ferroelectric materials due to 
their low, but finite mobility. In the orientational picture, oxygen vacancies, when adjoint 



2 



to an acceptor center, form electric and elastic defect dipoles in the bulk of a ferroelectric 

n 

domain [6| . The defect dipoles align parallel to the spontaneous polarization Pg by diffusion 
of the mobile oxygen vacancy in cage motion. Because of the relatively slow oxygen vacancy 
motion Q], the polarization directions of the aligned defect dipoles stay constant when the 
direction of Ps changes for short times. In this case, the defect dipoles generate an internal 
electrical field E* which stabilizes the domain pattern by increasing the force constant for 



the reversible displacement of the domain walls 



271]. This relaxation model has been well 



developed [28|. It bears two insufficiencies, though, the time dependence of aging is not 
reproduced and the absolute values of clamping pressures are low 2^. A second point of 
view about the role of oxygen vacancies in aging is the formation of ionic space charges 



30 1 which was originally proposed to explain space charge effects in BaTiOs single crystals 



3l| . Ionic space charges are well known for highly doped positive temperature coefficient 



resistors based on BaTiOq 



32( 1 . For aging the mobile charge carriers move to charged domain 



faces or grain boundaries and compensate polarization. This leads to an asymmetric charge 
distribution whereby a voltage offset arises yielding the known shift or deformation of the 
ferroelectric hysteresis 33| . The clamping pressure on domain walls generated by these space 
charges has not been treated mathematically for periodic domain structures. 

This paper describes quantitatively the formation of space charges in single domains 
of a periodic structure and shows the development of the defect distribution inside the 
domain. An estimate of bending and clamping pressures on domain walls and a comparison 



to the orientational picture 



14j are given. Electrostatic clamping of domain walls through 



the formation of space charges is calculated to be two orders of magnitude stronger than 
clamping through aligned defect dipoles for the same concentration of charge carriers. 

The model is independent of the type of point defect, as long as a diffusion constant can 
be assigned and the defect is charged. It can thus be equally well applied to hopping of 
electronic carriers. The oxygen vacancy was chosen for the numerical examples in order to 
be comparable to previous work, but does not preclude a statement on the physical nature 
of the mobile carrier. 
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II. GENERAL MODEL 



In order to study the effect of migration of charge carriers on aging, we chose a two- 
dimensional periodic array of domains cut by the grain surface, z = 0, perpendicular to 
the direction of spontaneous polarization which is along the z axis in Fig. [T] This model 
configuration is well-known in the physics of polarized media and was used for the study 



35l | and ferroelectric 
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37| 



of equilibrium and dynamic properties of ferromagnetic 
materials. We assume for simplicity an isotropic material of the grain occupying the area 
z > characterized by the relative dielectric constant Sf. The dielectric material outside 
the grain is assumed isotropic too and is characterized by the relative permittivity e^- As we 
previously showed by finite element simulation, the electric fields arising due to spontaneous 
polarization in a periodic multi-domain grain of finite dimensions generate a nearly perfect 
periodic pattern except for the very edges of the grain [29|. We thus consider the periodic 
domain array of Fig. [1] occupying the semi-space — oo < x < oo, z > as a representative 
model for a multi-domain grain of domain width a much smaller than the typical grain 
size. Due to polarization, the domain faces ai z = are alternatively charged with the 
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FIG. 1: Scheme of a 2D-array of 180°-domain walls crossing the grain boundary at a right angle. 
Straight arrows show the direction of the polarization and curved arrows the approximate directions 
of the local electric fields. 



bound surface charge density a = |Ps|, the spontaneous polarization value. If the length of 
the domains L along the z-axis is much larger than their width a along the x-axis, which is 
typically the case in experiment, field lines are effectively closed at the same side of the grain 
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(see Fig. [T]). As a consequence, both components of the electric field E°(x, z) induced by the 
alternating surface charge exponentially decrease towards the interior of the grain along the 
2;-axis 29|. Migration of the charged defects driven by the field E°(x, z) is then expected to 
occur in a volume near the grain surface. The domains may therefore be considered infinitely 
long along the ^-axis without introducing a substantial error. 

Let us now derive the equations of evolution for the charge and field distributions in 
the considered system. At any time t, the electric field E{x,z,t) is determined by the 
charged faces of the domains and the imbalanced charge density of free carriers p{x, z, t) = 
Qf [c{x, z, t) — Co] through Gauss' law 

VE = ^(c-co) (1) 

eoEf 

where c{x,z,t) is the local concentration of mobile carriers of charge g/, cq is the back- 
ground concentration of the immobile charge carriers of opposite polarity warranting total 
electroneutrality, and Eo is the permittivity of vacuum. In the initial state, the system is 
locally neutral assuming c{x,z,0) = Cq, the electric field Ei{x,z,0) = E°(x, z) is yet to be 
found. 

In presence of an electric field and a gradient of concentration, the fiow of charge carriers 
is given by the sum of drift and diffusion contributions to the particle current density: 

s = /icE - DVc (2) 

where and D are the mobility and diffusivity of charge carriers, respectively. We assume, 
for simplicity, that the latter two quantities are isotropic and connected by the Einstein 
relation, D = fikQ/qj with k the Boltzmann constant and G the absolute temperature. 
Migration of charge carriers is governed by the continuity equation: 

dtc = -V(/xcE) + DAc. (3) 

For boundary conditions to the system of equations ([1]) and ([3]) we assume chemical and 
electrical isolation of the grain. The first requirement means vanishing particle current 

Sz = ncE^ - DdzC = 0, (4) 

at the grain boundary, z = 0. The second requirement means vanishing total electric current, 

qfificE, - Dd,c) + eoSfdtE^ = 0, (5) 



at z = which results in a constant value of the ^;-component of the electric field at the 
grain boundary, dtEz{x,0,t) = 0. 

Eqs. ( HfS]) together with Eqs. (11151) complete the statement of the problem of charge 
segregation in a ferroelectric grain. In the next section we will observe how the system 
relaxes according to the equations of evolution (IT|51) . 

III. SOLUTION OF THE EQUATIONS OF EVOLUTION 

In this section we first calculate the field E''(x, z) in the virgin state of the system before 
the process of charge segregation starts. Then we formally solve equation ([1]) and find 
the total electric field E(x, z, t) for an arbitrary right-hand side. Finally, using the latter 
result, we numerically solve equation ([3]), self-consistently describing drift and diffusion of 
the mobile charge defects in the domain arrangement of Fig. [H 

A. Electric field in the virgin state of a multi-domain grain 

To use the bilateral symmetry of the problem, the origin is chosen in the center of the 
positively charged domain face as shown in Fig. [21 




FIG. 2: Scheme of expected charge redistribution induced by the local electric field within the 
central region —a < x < a, z > 0, presenting repeating element of the periodic domain arrangement 
of Fig. [TJ Thin layers of positive charge carriers piled up at the negatively charged domain faces as 
well as a wide (shaded) area depleted of mobile charge carriers near the positively charged domain 
face are shown. 

The bound charge density of the domain faces is represented by an alternating function 
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Pb{x,z) = a6{z) 



oo 
?i=— oo 



-1)"6' (^^-an + x^9(^^ + an 



X 



(6) 



where 5{z) and 6{x) are the Dirac 5-function and the Heaviside unit step function, respec- 
tively. The electrostatic potential induced by this bound charge is given by the expression: 



1 



2TTeo{ef + ed) J- 



dxo / dzoPb{xo,zo)lia 



X — Xo 



a 



+ 



Z- Zo 



(7) 



in both areas 2; > and z < 0. The formula ([7]) is simply a superposition of the potentials 
generated by straight parallel charged lines located at the grain boundary z = between 
the media with dielectric constants Ed and Ef [35|. 

The 2;-component of the electric field created by the bound charge, E*^ = — Vy^^,, may be 
directly calculated by substitution of pb, Eq. ([6]), into Eq. ([71) and subsequent summation 
which results in the form 



E'^Jx, z) 



2a 



■ arctan 



cos(7rx/a) 
sinh(7r2;/a) 



(8) 



7TEo[Ef + Ed 

valid for both media. 

Direct calculation of the other field component, = —dx^Pb, is more complicated because 
of slow convergence of the appropriate series. Instead, E^ may be calculated for 2; > 
from Gauss' law VE*^ = 0, taking into account that, from the symmetry of the problem, 
E^{0,z) = E^{±a,z) = 0. Proceeding with integration of the latter Gauss' equation over 
distance along the x-axis and using the mentioned boundary conditions one finds the form 

cosh(7r2;/a) + sin(7rx/a 



E^Jx^z) 



a 



7rEo{Ef+Ed) 



In 



cosh(7r2;/a) — sin(7ra;/a) 



(9) 



valid inside and outside the grain. Both field components exhibit periodic dependence along 
the expected from the periodic domain arrangement, and exponentially decrease 

at large distance from the charged surface \z\ ^ a, as expected from the previous finite 



element simulations 



29| . We note here that the short range fields, Eq. f lHllQl) . may be very 



arge. For example, in BaTiOa the depolarization field amplitude may be as high as 10^ V/m 



291] . The presence of high local fields were confirmed at least partly in observations on the 



acceptor doped BaTiOa, where internal fields were experimentally found in the range of 
10^ V/m for up to 1 mol% Ni-doping and temperatures below 80°C jc]. 
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B. Electric field due to redistribution of charged defects 



At any arbitrary moment, the total electric field in the system may be conveniently 
decomposed as E = E'^ + E*, where the field E" is determined by the charged faces of 
the domains, Eqs. and the field E' is generated by the distribution of the charge 

density in the area z > 0. Thanks to periodicity and the bilateral symmetry of the initial 
conditions, both the charge density and the electrostatic potential remain periodic and 
bilaterally symmetric in the course of the charge redistribution, as illustrated in Fig. [2l This 
allows us to consider the region —a < x < a as a repetitive basic unit of the system and 
confine ourselves to the consideration of processes in this area. To get a full description of 
the electric field under these circumstances, it is sufficient to construct the Green's function 
of the symmetrical Neumann problem in the mentioned region, Gs{x, z\xo, zo), so that the 
electrostatic potential induced by redistribution of charged defects with Zq > may then be 
presented in a form {s^ 



ipi{x,z,t) 



dxo / dzop{xo,zo,t)Gs{x,z\xo,zo), 
Jo 



(10) 



followed by the field expression E* = —Vipi. 
The Green's function satisfies the equation 

1 



/\Gs{x,z\xq,Zo) 



EoEf 



5{z - Zq) [S{x - Xq) + 5{X + Xq)] 



(11) 



with boundary conditions d^Ggix = ±a, z\xq, Zq) = 0. The latter requirement is a conse- 
quence of the constraint Ex{±a, z) = inherent to the chosen domain arrangement. Bound- 
ary conditions for the electrostatic potential, Eq. ( JTOl) . on the interface between the two 



media at ^ = 



39j impose two additional boundary conditions on the Green's function 



Gs{x, -0\xo, zo) = Gs{x, +0\xo, zq) 

Edd^Gsix, -0\xo, Zq) = Efd^Gs{x, +0\xo, Zq) 



(12) 



Using the fundamental solution of the 2D Poisson equation 



39| and taking into account 



periodicity of the problem the solution of Eq. (ITTI) may be reduced to summation of a series: 



Gs{x, z\Xo, Zq) 



2TTEo{Ef + Ed) 



E 



In 



X — Xq 



a 



2n 



-Xq) 
(13) 
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for the area 2; < and 



1 



Sf + Sd 



4:71606 f 

X — Xq 



In 



X — Xo 



2n 



a 



2n + 



Z + Zq 



+ {xo 



-Xo) 



(14) 



for the area z > 0. 

Because of slow convergence of this series it is more convenient to perform summation 
for the derivatives dxGg and dzGg and then to restore the function Gs itself by integration 
using boundary conditions. This leads after all to the solution of Eg. (11 II) 

1 



Gs{x,z\Xq,Zq) 



2-K6Q{6f + 6d) 

for the area z < and 



In 



7r(2-2o) ^{x-Xq) 

cosh COS 



+ (xo -Xo) (15) 



Gs{x,z\xq,Zq) 



1 



In 



cosh 



. In 

ATX6o6f 

TC{Z + Zo) 



cosh 



tt{z - Zo) 



COS ■ 



7r(x — Xo) 



cos ■ 



vr(x - Xo) 



(xo 



-XoJ 



(16) 



for the area z > 0, which are periodic, bilaterally symmetric and satisfy the proper bound- 
ary conditions. Now, from the expressions fll0|15lll6p . the components of the electric field 
induced by the redistribution of charged defects may be obtained. It is easy to verify that 
the total electric field satisfies the boundary condition Ex{x = ±a, z) = for any bilaterally 
symmetric charge density p{x,z,t). 



C. Numerical solution of the evolution equations 

Having solved equation ([1]) explicitly allows for the implementation of a simple direct 
Euler scheme for numerical treatment of the problem. Space and time will be discretized. At 
every time step, the change in the carrier concentration will be calculated from the previous 
values of the concentration and the electric fields using Eq. ([3]). Then, the updated values 
of the field will be calculated directly from the updated values of the concentration using 
Eq. ffTOj) . The calculation is repeated until convergence. Taking into account the bilateral 
symmetry of the problem, it is sufficient to consider the charge redistribution within the 
area < x < a. 
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We first introduce dimensionless variables which is helpful for the following numerical 
analysis and reveals those parameters of the system which are relevant to the relaxation 
process. Dimensionless coordinates are naturally introduced as X = x/a and Z = z/a. The 
dimensionless field F = E/i?* is expressed in units of the characteristic value E* = a/26o6f. 
The system reveals two characteristic time scales: the drift time = a/ fiE* and the diffusion 
time td = a? /D. For the typical parameters involved, td ^ therefore we will introduce 
dimensionless time as T = t/r^. The concentration of defects is now reduced to n{X, Z, T) = 
c{x,z,t)/c* with the characteristic value c* = a/2aqf. The latter has the physical meaning 
of a concentration of defects on an area a^, which completely neutralizes the bound charge 
a at the domain faces. The reduced initial concentration uq = Cq/c* measures whether the 
density of defects is high or low with respect to the charge compensation concentration. 

The continuity equation ([3]) now acquires the form 



where all differentiations are performed with respect to the dimensionless variables. The 
parameter (3 = t^/tij ^ 1 characterizes a weak contribution of diffusion to the migration of 
defects in ferroelectrics. It is now seen from Eq. (1171) that only two composed parameters, 
hq and (3, control the relaxation process. 

Though the parameter /3 may be rather small, it cannot be neglected as is clearly seen 
from the boundary condition for the particle current, Eq. (jl]), taken in a dimensionless form 



otherwise this boundary condition is not compatible with the initial conditions. The finite 
value of P means compensation of the drift contribution to the current by the diffusion 
contribution at the grain boundary and this way defines the structure of a thin layer of 
charged defects piling up at this boundary. 

Eq. ( |T7I) is supplemented by expressions for the dimensionless field F = F*^ + F* which 
can be easily derived from Eqs. ( |8|9l) and ( fT0]|15|16l) . namely. 



dtu = —n{n — no) — FVn + f3An 



(17) 



nFy - (3dyn = 0, Z = 0, 



(18) 



F'AX,Z) 



1 2ef , rcoshvrZ + sinvrX 

■' 1 v-i 



iref + Sd cosh ttZ — sin ttX 




{X,Z) 



2 2ef [costtX 

— arctan — - — — 

TT Ef + Ed [smh ttZ 



(19) 
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FUX, Z,T)= / dXo / dZo Z\Xo, Zo) [n{Xo, Zq, T) - n,] 



and 

f;,(x,z,t)= ) 



where the kernels in this integral are presented by the functions 

Ef sin7r(X — Xq) 



(20) 



fx{X, Z\Xq, Zq) 

fz{X, Z\Xo, Zq) 



2{ef + Ed) cosh7r(Z — Zq) — cos7r(X — Xq] 
Sf sinh 7r(Z — Zq) 



cos7r(X — Xc 



2{£f + £cl) cosh7r(Z — Zq) 
for Z < and by functions 

Ef-Ed sin7r(X-Xo) 



+ (Xo- 



-Xo), 

. -Xo) (21) 



fx{X, Z\Xq, Zq 



Ef + Ed cosh7r(Z + Zq) — COS7r(X — Xq 
sin7r(X — Xq) 



cosh7r(Z — Zq) — cos7r(X — Xq 



+ iXo 



-Xn 



fz{X, Z\Xq, Zq 



+ 



Ef - Ed 



sinh 7r(Z + Zq 



Ef + Ed cosh7r(Z + Zq) — cos7r(X — Xo^ 
sinh 7r{Z — Zq) 



+ {Xo 



-Xq 



(22) 



cosh7r(Z — Zq) — cos7r(X — Xq) 
for Z >0. 

Since the system remains electrically neutral within the domain of integration during the 
redistribution of defects, arbitrary constants may be added to the kernels fl21|22l) without 
changing the results of integration in Eqs. fl2UI) . This property is used in the numerical 
procedure to facilitate the conversion of the integrals in Eqs. (!20l) . 

As an example, we now consider the aging process in BaTiOa. For the numerical simu- 
lations, the material parameters of BaTiOs at room temperature are taken from Wernicke 
and Jaffe et al. |4o|, kli, namely, = 2.71 ■ 10"^ C/cm^, Ef = 170, fx = 1.73 ■ lO^^o mVVs, 
a = 0.5/im and gj twice the elementary charge, implying positively charged oxygen vacancies 
as mobile defects. For the dielectric medium between ferroelectric grains we take the same 



but non-polarized material with Ed = 170. This yields c* = 1.69 ■ 10 cm 



1.61 -lO^ s. 



To = 1.14 ■ 10^ s. As was shown in one dimensional simulations 29(], the parameter (3 << 1 
has no effect on the dynamics of the relaxation. The only physical characteristic depending 
on /3 << 1 is the thickness of the positively charged layer of defects piling up at the negative 
face of the domain. To make this layer visible in figures and to avoid numerical problems 
invoked by the strong gradients of the defect density we take the value /5 = 5 ■ 10~^ instead 
of the actual ratio t^/td = 1.4 ■ 10~^ for our simulations. 
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FIG. 3: Distribution of oxygen vacancies cy^(X, Z) over the reference area 0<X<1,0<Z<4 
at time T = 0.05 for an initial concentration of defects cq = c* = 1.69 • 10^^ cm~^. 

A snapshot of the development of the defect concentration profile over the reference area 
0<X<1,0<Z<A starting with the background defect concentration no = 1 is presented 
in Fig. [3] for the moment T = 0.05. A wide depleted zone forms near the positively charged 
face at < X < 0.5, Z = and a very thin excess charge layer of high concentration near 
the negatively charged face at 0.5 < X < 1, Z = 0. 
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FIG. 4: Defect concentration profile along the line X = 0.25 for a succession of times T = 
0.1, 1, 2, 3, 4, 5 (from left to right) for the initial concentration of defects cq = no • c*. 
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FIG. 5: Defect concentration profile along the line X = 0.75 for a succession of times T = 
0.1, 1, 2, 3, 4, 5 (from left to right) for an initial concentration of defects cq = uq ■ c*. 

The structural difference of these two space charge areas is better seen in Figs. H] and 
|5] presenting vertical cross sections of the concentration profile along the lines X = 0.25 
and X = 0.75, respectively. A succession of snapshots of the concentrations along the line 
X = 0.25 (Fig. H]) exhibits the evolution of the charge defect density near the positive face of 
the domain. The profile positions at the moments T = 4 and 5 cannot be discerned any more 
which indicates saturation at time T ^ 5 (corresponding to t ~ 8 ■ 10^ s). The characteristic 
width of this zone in the final relaxed state is of the order of unity. The defects piling up 
near the negative face of the domain form a much thinner layer of a characteristic width of 
the order of f3 as is represented by concentration profiles along the line X = 0.75 in Fig. O 
The final relaxed state is reached also at about T ~ 5. The corresponding evolution of the 
front cross section of the concentration profile along the line Z = shown in Fig. exhibits 
saturation at about T ~ 5, too. 

In our model, drift-dominated migration of the charged defects is caused by local electric 
fields near the charged faces of a grain. This migration process only stops, if either no mobile 
defects remain in the area where fields are present or there is no remaining field in the area 
where the defect concentration is not zero. The process of field compensation due to defect 
migration is exemplified by the evolution of the electric field component = + at 
the line Z = 2 represented in Fig. [3 F^ saturates at the values opposite to the local values 
of the initial electric field F^ determined by the bound surface charge. Relaxation leads to 
an energy minimum where the system will resist any change of the domain wall positions. 
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FIG. 6: Defect concentration profile along the line Z = for a succession of times T = 
0.1, 1, 2, 3, 4, 5 (upwards) for the initial concentration of defects cq = uq ■ c*. 
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FIG. 7: The electric field component plotted along the line Z = 2 for a succession of times 
T = 0, 0.1, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 in a system with the initial concentration of defects cq = c*. 
The arrows show the direction of evolution. 



The final distribution of free charges then determines the equilibrium domain configuration 
of the system. For an effectively low mobility of the free charge carriers the transition above 
the Curie point will not readily rearrange the charge carrier configuration due to thermal 
excitation. The defect charge density then determines the subsequent domain configuration 
after re-cooling the sample to low temperature. Experimentally it is observed that the 
original domain configuration is reproduced to a great extent [3, a, l42 |. 
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IV. FORCES EXERTED UPON A DOMAIN WALL 



From the known development of the charge density and the electric field in our model, 
the time dependent forces exerted upon domain walls can be evaluated. Using the general 
formula derived by Nechaev et al. [431] and taking into account only electrostatic contribu- 
tions to the energy one can obtain the local pressure / exerted upon a wall. For a straight 
rigid wall considered here, one finds / = 2PE where P and E are the local values of sponta- 
neous polarization and electric field, respectively. This relation is reduced, in the geometry 
of Fig. [H to / = 2PsEz (note that, in the case of the same arrangement of the 90°-domain 
walls, the a would merely decrease by a factor of -\/2 and the force by a factor of 2, the 
configuration and results are otherwise identical). 

Although only one end of the domains is present in the mentioned geometry of Fig. [1] 
it is obvious that similar segregation of the charged defects occurs on the other end of the 
domain, too. This results in the antisymmetric force of opposite sign exerted upon the 
domain wall on the other end of the domain yielding a total force equal to zero. This force 
cannot move the domain wall as a whole or prevent its motion but it may lead to bending of 
the wall violating our assumption of rigid straight domains. This is frequently encountered 
in real systems. Domains forming needle tips near external interfaces are commonly observed 



44 



45|. In this case, part of the compensation arises within the bulk and not only right at 



the grain interface. The final defect distributions will be different from the case calculated 
here, but the essential effect of bending will remain the same. Our model of drift of free 
charge carriers also supports a coalescence of domains rather than their splitting. Without 
any further details included in the model, it contradicts the experimental observations of 
Ikegami ad Ueda of domain splitting during aging [lo| . 

The evolution of the bending pressure f{T) averaged over half the domain wall length, 
assumed as long as L = 20a, is shown in Fig. [8] for three different values of the initial 
background concentration of defects. It is seen that systems with smaller concentrations 
need an inversely longer time to relax. For the system with no = 1 it takes about T ~ 5 
while for the system with no = 0.5 this time is roughly doubled. All curves can be well fitted 
by the exponential form /otanh (anoT/2) where the parameters /o — 1 MPa and a ~ 1 and 
slowly increase when no decreases. A reliable simulation of defect concentrations smaller 
than uq = 0.5 is impossible on the chosen template (0 < X < 1, < Z < 4) since in this 
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case migration involves defects from a wider area in order to compensate the bound charge 
at the domain faces. 




FIG. 8: Bending pressure / as a function of time T = t/r^ at room temperature for BaTiOs is 
plotted for three different sample concentrations of oxygen vacancies cq = uq ■ c* with uq = 0.5, 1 
and 2 (solid lines). Dashed lines show fitting of the pressure by the function /o tanh {anQT/2) with 
parameters /o = 1.095, 0.91, 0.66 MPa and a = 0.91, 0.86, 0.76 for no = 0.5, 1 and 2, respectively. 

One more general feature of time dependencies of the bending pressure in Fig. [8] is worth 
discussion. All the curves exhibit a small region at small times where the value of pressure is 
negative. This is not an artefact of the numerical discretization procedure but has a physical 
meaning. Indeed at any time, the characteristic width of the positive space charge zone near 
the domain face 0.5 < X < 1, Z = is of the order of f3 which follows from the boundary 
condition, Eq. (jlj). On the other hand, at the very beginning of charge defect migration, the 
characteristic width of the negative space charge zone in the area < X < 0.5, Z > is less 
than f3. This means that a negative value of the field component Fz prevails at the domain 
wall at this stage. This is confirmed by the dependence of Fz on position Z for different 
times as presented in Fig. M 

The above considered bending force does not directly describe the aging phenomenon 
as long as rigid straight domain walls are retained. In fact, the total force exerted upon 
the walls remains equal to zero during the defect redistribution if both ends of the domains 
are taken into account. Nevertheless, the loss of domain wall mobility characteristic of 
aging may be captured in this model, too. Indeed, the segregation of charge carriers in the 
fixed domain framework of Fig. [1] leads to the relaxation of the energy of the electrostatic 
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FIG. 9: Snap-shots of the distribution profile of the field component along the domain wall 
for the succession of times T = 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1 (upwards) for BaTiOg at 
room temperature for an oxygen vacancy concentration cq = uq- c* . 

depolarization field. The decrease of this energy per unit length of domain wall measures 
the clamping pressure preventing the displacement of the wall from the energy minimum: 



The dependence of this pressure along the length of the wall is shown in Fig. [10] for a 
succession of times. The magnitude of the pressure saturates as expected at about a time 
T ~ 5 for a defect concentration of Cq = c*. The corresponding peak value of the pressure 
around 1.5 MPa is comparable with the average bending pressure at the wall, Fig. [HI The 
magnitude of the saturated pressure increases monotonously with the defect concentration 
Co as is seen from Fig. [HI 

The irreversible migration of charged defects entails growing immobilization of the domain 
walls and, consequently, enhancement of the coercive field, Ef.. To estimate this effect one 
should compare the pressure ~ PgS exerted by the external field, upon a domain wall 
with the clamping pressure, Eq. (!23l) . averaged over the domain wall length, L. This results 
in the following estimate for the coercive field 



J^sJ^ Jo 

where integration over the half-length of the wall accounts for the other end of the domain. 
Evaluation of the time-dependent coercive field assuming the typical length of the domain 




(23) 




■L/2 



(24) 
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FIG. 10: Snap-shots of the clamping pressure distribution along the domain wall for the succession 
of times T = 0, 0.1, 1, 2, 3, 4, 5 (upwards) for BaTiOa at room temperature for an oxygen vacancy 
concentration cq = uq ■ c*. 
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FIG. 11: The saturated clamping pressure distribution along the domain wall for BaTiOs at room 
temperature for the oxygen vacancy concentrations cq = no • c* with no = 0.5, 1, 1.5 and no = 2 
(upwards) . 



wall L = 20a obtains a characteristic value of — 1 kV/cm which is of the order of 
the coercive field in unaged bulk samples of BaTiOa In fact, the magnitude of the 

clamping pressure and, consequently, the value of the coercive field may be substantially 
larger then it was estimated using Eq. fl2^ . Firstly, the peak value of the pressure has to be 
approximately doubled if one takes into account the reduction of the energy of electrostatic 
field in the dielectric material outside the grain which is approximately the same as in the 
ferroelectric area assuming ea = Sf. Secondly, the consideration of the anisotropy of the 
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dielectric constant is expected to scale up the pressure together with the energy gain by the 
factor of \/ Eal^c which is about 6 for BaTiOs. Finally, values of few MPa are expected 
for the average clamping pressure at the domain wall and the values of few kV/mm are 
expected for the coercive field due to charge carrier migration which is in agreement with 

on 

the characteristic values observed on the aged samples of BaTiOs p, |8|. Accordingly, the 
coercive field, Eq. (p^ . multiplied by the factor of 12 is shown in Fig. [12] in physical units 
to compare with known experimental data. The dashed line shows that the time behavior 
of Ef. mimics logarithmic time dependence for durations less then a few r^. 




FIG. 12: Coercive field due to charged defect migration as a function of time (solid line) for the 
oxygen vacancy concentration cq = • c*. Dashed line shows fitting with logarithmic dependence 
for intermediate times. 

One more essential factor which brings about enhancement of the coercive field is that 
the minimum energy of the system will further substantially decrease if domain wall bending 
is allowed contributing to the increase of the clamping pressure, Eq. ( l23l) . This mechanism 
is, however, beyond the consideration in our model of rigid walls. 

The above obtained values are much larger than typical magnitudes of clamping pressure 
arising due to dipole re-orientation jl^. Indeed, for uniformly aligned dipoles in the latter 
mechanism, the dipole moments exert upon the domain wall a clamping pressure Por — c^E^d 
where the dipole moment d = qfl/2 with the dipole length of / = 4- 10~^''m. For the material 
parameters assumed in the above estimations and cq = c* this results in the peak value of 
the clamping pressure Por = 9.7 kPa which is two orders of magnitude smaller than that in 
the drift mechanism. 
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A common feature of these two aging mechanisms is that both dipole re-orientation and 
defect migration occur in those areas where a depolarization electric field is present. In 
respect thereof these mechanisms can be classified neither as volume nor as boundary ones 
as was suggested in the recent works by Zhang and Ren [3, O] but rather as geometry 
dependent. Indeed, in the two-dimensional periodic array of domains considered here the 
depolarization field is present only near the grain boundaries causing both defect drift and 
dipole re-orientation only in this area. On the other hand, in the single domain state of a 
Mn-doped BaTiOs single crystal observed in Ref. Sj a one- dimensional geometry is virtually 
realized where the depolarization field is present in the whole sample 291] and invokes both 
dipole re-orientation and defect drift in the whole volume. 

V. CONCLUSIONS 

In this work, we have considered migration of charged defects as a possible reason for 
aging in ferroelectrics. The model is based on two main assumptions: 1) existence of mobile 
carriers of ionic or electronic nature in the bulk material and 2) presence of strong local 
depolarization fields due to bound charges at the domain faces. The first assumption is 



based on direct measurements of the conductivity in perovskites 26|, the second one was 
at least partly confirmed in observations of Ref. j^. Solving self-consistently the drift- 
diffusion equation together with the Gauss equation for the fixed two-dimensional domain 



array 
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35l . l36l . |37] reveals gradual formation of space charge zones compensating the field 
generated by charged domain faces. Charged domain walls, which are tip-to-tip or tail-to- 
tail configurations of the polarization in adjacent domains, are electrically equivalent to our 
model. The biggest difference arises due to the fact that charged domain walls are often 
observed as needle tip domains in single crystals [45]. The geometry is thus considerably 
different from the model of parallel domain walls presented in our paper, where only periodic 
straight domain configurations are captured. 

The process of charge defect migration is accompanied by the reduction of the energy of 
the electrostatic depolarization field which leads to the energy minimum where the system 
will resist any change of the domain pattern. The characteristic time of this relaxation 
depends on the doping and is typically about 5-r^~8-10^s~9 days, where is a time 
of drift over the distance of domain width. That is why, after aging, a clamping force at 
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a domain wall arises if an external electric field attempts to shift the domain wall from its 
initial position. This force may be estimated from the calculated energy gain due to the 
reduction of the depolarization field. The peak value of the clamping pressure is in the 
range of 1 ^ lOMPa but the pressure is distributed very inhomogeneously along the domain 
wall concentrating near the domain ends. Nevertheless, the total value of the clamping 
force at the domain wall results in the characteristic coercive field of few kV/mm which is 

n 

comparable with that observed on the aged samples of Mn-doped BaTiOa [8|. 

Clamping pressures on domain walls in the presented two-dimensional model are consid- 
erably lower than in the uniaxial case [29:| and approach macroscopically observable values. 
They are two orders of magnitude larger than in the picture of defect dipole re-orientation 
fl^ and are thus a plausible mechanism for aging in ferroelectrics. In contrast to the one- 
dimensional case with only one characteristic value of electric field, Ed = Ps/sfSo, treated 
earlier [2^ the two-dimensional model exhibits seemingly a wide spectrum of characteristic 
times according to the position-dependent values of the electric field E(x,?/). This allows 
one to expect a time dependence of the clamping pressure in a two-dimensional array of 
domains different from the one-dimensional case (29i| . Nevertheless, comparing time evolu- 
tion of the field and defect concentration in Ref. 29|] with Figs. pi[5]|6|7|) one observes a 
striking similarity between them. We are thus concerned with a single characteristic time 
constant r^. = t^/uq characterizing the relaxation of the system. This time is independent 
of the width of the domains, a. In fact, = EfEo/X with A = g/Co/i being the conductivity 
of the material is the Maxwell relaxation time which only depends on the mobility and local 
concentration of the mobile carriers. This in turn means that a distribution of grain sizes in 
the material and accordingly a distribution of domain sizes does not entail a distribution of 
characteristic relaxation times. The logarithmic time dependence of the dielectric constant 
during aging yet remains to be explained. 

A crucial parameter for the plausibility of the time scale in our simulations is the mobility 
of charged species in a ferroelectric material. The mobility of oxygen vacancies considered is 
still a highly disputable issue. The activation barrier for this ionic defect is usually estimated 



3,Q, 



in the range of 0.9-1.1 eV in both experimental works and first principle calculations 

which makes the migration of oxygen vacancies over the distance of the order of the 



domain width ~ 0.5 um most unlikely. On the other hand, the estimations o: 



in the Refs. 26, 



30, 



50, 



5l| are similar to or higher than that given in Refs. 



the mobility 



4l| which 
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we used for simulations in our study. We would like to stress here therefore that the nature 
of the charge carriers plays no important role for the model presented. These may be also 



electronic carriers as was suggested in Refs. 47|, |52] • In any case our input parameters agree 



with direct measurements of the conductivity of perovskites indifferent to the nature of the 



charge carriers 26 1. 



It is evident that any real system will contain more than one mobile charge carrier. In 
case their mobilities or concentrations are considerably different, the final distribution of 
defects of the more mobile / more frequent carrier will determine the field environment 



for the drift of the second carrier as it was discussed in Ref. [47l]. The solutions from the 
present calculation would have to be taken as starting condition and iteratively the final 
solution could be found. In case of equal mobilities and concentrations, a coupled system 
of equations has to be solved which is the issue of forthcoming work. Similarly the local 
potential wells for the domain wall which determine the dielectric constant will be given in 
a future publication. 
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